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We study the three-dimensional Ising model at the critical point in the fixed-magnetization 
ensemble, by means of the recently developed geometric cluster Monte Carlo algorithm. We define 
a magnetic-field-like quantity in terms of microscopic spin-up and spin-down probabilities in a given 
configuration of neighbors. In the thermodynamic limit, the relation between this field and the 
magnetization reduces to the canonical relation M(h). However, for finite systems, the relation is 
different. We establish a close connection between this relation and the probability distribution of 
the magnetization of a finite-size system in the canonical ensemble. 
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d . I. INTRODUCTION 

C/3 . 

Second order phase transitions display many interesting and subtle properties associated with scale invariance and 

£H universality at critical points. Some of these, such as power-law singularities of the free energy and other quantities 

at the critical point, and their critical exponents and amplitudes, have been studied rather thoroughly (see, e. g., jl|). 
Among less investigated items are the universal characteristics of finite-size effects. These are important for analysis 
of experiments with finite samples, as well as for computer simulations, which necessarily have to deal with finite-size 

O ' systems. 

O , In clear distinction from, for example, critical indices, the finite size effects depend crucially on the nature of the 

statistical ensemble under consideration. To be concrete, let us consider one of the standard model systems of the 
phase transition theory — the three-dimensional (3D) Ising model on the simple cubic lattice, with nearest-neighbor 

^ • interactions. 

According to universality, this model describes the critical properties of a wide range of systems of a different physical 
nature, including second order phase transitions in uniaxial magnetic systems and the liquid-vapor critical point. The 
statistical ensemble most commonly used (usually denoted as canonical ensemble) is defined by the partition function 

OS', Z c (h) = y^expj/gy^s^j + fty^Stj, Si = ±l, (1) 

O^ . {s,} (ij) i 

-f— > I 

where the sum includes all the 2 W possible configurations of a total number of N spins. This ensemble is perfectly 

^ • natural for applications to magnetic phase transitions, where s, = ±1 corresponds to physical spin at the lattice site 
i pointing up or down, respectively. However, in the language of the lattice gas, Si = ±1 corresponds to an occupied 
or unoccupied site, respectively: the total number of particles fluctuates. The canonical partition sum Eq. fll|) in 
the spin language thus corresponds with the grand partition sum in the lattice gas language. To avoid confusion, we 

CJ i emphasize that in the rest of this paper we are using the word 'canonical' in the language of the Ising model: it means 
that the magnetization is allowed to fluctuate freely. 

In contrast, many real experiments and simulations of liquid- vapor systems are performed with a fixed number of 
particles. Fixing the number of lattice-gas particles is equivalent, in the language of magnetic systems, to fixing the 
total magnetization M to tai = Sj=i s « °f the system, or, in other words, fixing the average magnetization per spin 



Thus we will be interested in the properties of the 3D Ising model in the fixed-magnetization ensemble, 
Note that the magnetic field h is absent; it can only contribute a constant factor. 



One of the main difficulties encountered in computer simulation studies of critical phenomena is the critical-slowing- 
down phenomenon. For a number of spin models in the canonical ensemble this problem has been largely overcome 
with the invention of cluster algorithms [g,|| . Until recently, no similarly useful algorithm has become available for the 
fixed-magnetization ensemble. This situation has now changed by the development of a geometric cluster algorithm 
P,|l . We have used this algorithm extensively in this work to efficiently simulate systems of fixed magnetization at 
the critical point. 

In the canonical Ising model described by Eq. (fy), the magnetic field h is an adjustable external parameter. In 
contrast, the magnetization is a fluctuating observable. For each configuration taken from the ensemble one can sample 
its magnetization per spin M = -A> X)i=i s i- Having accumulated M over a sufficiently large set of configurations, one 
can construct the probability distribution P(M) p|-[Tc|], and determine various expectation values such as (M), {M 2 }. 
Examples of such probability distributions at the critical point are shown in Fig. EL 

On the other hand, for systems in the fixed-magnetization ensemble described by Eq. (g) , the roles of h and M are 
interchanged: now M is the adjustable parameter, and it is intuitively clear that there should be some way to define 
an observable, which we denote by h (to avoid confusion with h in Eq. (|lj)), that will correspond to the magnetic field. 
Thus h will be a fluctuating quantity, that can be sampled on a microscopic level from configurations taken from the 
fixed-magnetization ensemble. In the limit N — > oo in both ensembles (such that the correlation length vanishes in 
comparison with the system size), the fluctuations in M and h become negligible, and the difference between h and 
h vanishes. 



In Section O we discuss the definition of h and its properties. In Section III we establish the relation between 
the function h[M) in the fixed-M ensemble and the probability distribution P\M) in the canonical ensemble. We 
conclude with a discussion of the relation between the function h(M) in the fixed-M ensemble and the function h(M) 
in the canonical ensemble, and with a summary of our main results. 

II. THE MAGNETIC OBSERVABLE H(M) FOR THE FIXED-A/ ENSEMBLE 

We will now describe a definition of h(M) that is based on statistical analysis of the local environment of a given 
spin. By local environment we mean the set of neighbors with which this spin interacts. In our particular case of 
the Ising model with nearest-neighbor interactions on the simple cubic lattice, the local environment consists of 6 
spins on the neighboring sites. The local environment has 2 6 possible configurations which divide in 7 types: 0+6— 
(zero spins up, 6 down), 1+5—, 2+4—, 3+3—, 4+2—, 5+1 — , 6+0— (six spins up, zero down). The simplest way to 
Monte Carlo sample h(M) is on the basis of the symmetric case: 3+3— ||. For every Monte Carlo configuration, go 
through all sites, and select all spins with the required 3+3— local environment. Then compute the average (sq) of 
the selected spins, and define h by 

(so) = tanh/i. (3) 

It is also possible to employ, instead of 3+3—, any other of the 7 types of local environment. In these non-symmetric 
cases the definition reads 

6 

(s ) = tanh(h + /?^ s 0i J , (4) 

»=i 

where the soi are the nearest neighbors of sq. 

One easily notices that the definition is constructed in such a way that h corresponds, on the mean field level, to 
the external field h in Eq. (En. 

Now it is interesting to see what are the results of Monte Carlo simulations for h(M). As has already been 
demonstrated in pLpl, a ^ the critical temperature h(M) practically coincides with the relation h(M) in the canonical 
ensemble as obtained by Monte Carlo simulations, provided M is sufficiently large, so that the correlation length is 
sufficiently small in comparison with the system size and the finite-size effects are suppressed. At the same time, the 
striking feature of h(M) for not-so-large M is its nonmonotonic behavior. First h(M) goes negative at small M, then 
it begins to grow, and finally assumes the usual behavior at larger M ||. This is clearly seen in Fig. J2J (diamonds), 
which shows Monte Carlo results obtained by means of the the geometric cluster algorithm ElO . 

In the remaining part of the paper, we will give the explanation of this behavior (which turns out to be a peculiar 
kind of finite size effect), by establishing a close relation between h{M) in the fixed-M ensemble, and the probability 
distribution P(M) in the canonical ensemble. 



III. CONNECTION BETWEEN H(M) IN THE FIXED-A/ ENSEMBLE AND THE PROBABILITY 
DISTRIBUTION OF M IN THE CANONICAL ENSEMBLE 

Considering the fixed- M ensemble, Eq. (g), one notices that it can be obtained by taking the canonical ensemble 
(M, and cutting from it the subset satisfying the constraint, ^ Si — NM . Within this subset we still have the usual 
Boltzmann probabilities exp{/3J2u j\ s i s j} f° r individual configurations. 

This makes it possible to establish a relation between h(M) in the fixed-M ensemble, and the properties of the 
system in the canonical ensemble. The definition of h(M) described in Sect. || is equivalent to the following. Let us 
take the fixed-M ensemble and concentrate our attention on one particular lattice site, and on the spin located there. 
Let us perform the following measurement. For every configuration consider the local environment of our selected 
site. If it is not 3+3—, do not measure anything for this configuration. If it is 3+3—, measure the selected spin Sq 
and store it. Finally, find (so), and use Eq. (^) to determine h. 

One notices that, as long as we are performing a thought experiment, we need not care about the Monte Carlo 
statistics. We can just stick to one site and get the same (so) without averaging over all sites, because they are 
equivalent. 

Up to now we have distinguished between 7 types of local environments, such as 3+3—. Let us go a bit further and 
treat separately all 2 6 possible local environments. In other words, the measurement of (sq) is now performed in an 
even smaller subset of the fixed-M ensemble: also the six spins forming the local environment of sq are fixed. In the 
case that the predetermined local environment is of the 3+3— type, it may seem that there is no interaction between 
So and the remaining system of N — 7 spins. Nevertheless, the fixed-M ensemble probabilities that so is +1 or —1, 
which we denote, respectively, by P+ and P_, are not equal in general. These probabilities may still depend on the 
magnetization of the remaining system, which is coupled to so by the overall magnetization constraint: 

6 N 

so+^2 soi + J2 Si = Yj Si - NM - ( 5 ) 

i=l ieRS i=l 

The total magnetization of the system is thus expressed as the sum of three terms: the local spin so, the sum of its six 
neighbors and the magnetization of the remaining TV — 7 spins denoted as X^eijs s i; where RS stand for "remaining 
system" . 

The conditional probabilities P± can be more explicitly written as 

P ± = P(s = ±l\NM,s i---s 06 ) (6) 

The two conditional arguments specify the total magnetization NM and the states of the 6 neighbor spins. We 
now make the connection with the canonical probabilities P c which include the magnetization as an unconditional 
argument. We use the zero-field canonical probabilities i.e., h = in Eq. (|l|). 

P± = P-^JVM |*oi ' ' ' s 06 )Pc(s = ±1, NM\s i • • • s 06 ) (7) 

We may slightly rewrite this by substitution of the probability P c by P c which is equal but uses the magnetization of 
the N — 7 remaining spins as its second argument: 

6 

P± = P-\NM\s m ■ ■ ■ s 06 )Pc(so = ±1, ^ * = NM ~ J2 s °* - s °l So1 ' ' ' s ° 6 ) ( 8 ) 

t£RS i=\ 

Let us first consider the simplest case 53j=i s oi — 0. Thus the canonical probability P c does not depend on its first 
argument, which can thus be skipped: 

P± = ip- x (JVAf |*oi • ' ' so 6 )P c ( J2 Si = NM T 1 ' s ° 1 ' ' ' s ° 6 ) ( 9 ) 

Therefore, 



i£RS 



P+ = P c{E l& RS s l = NM- 1|SQ1 ■ ■ ■ s 06 ) 
P - ~ PciEieRS s l = NM + l|s i ■ • • s 06 ) ' 



(10) 



The condition soi • ■ • so6 in effect introduces a defect in the remaining system: an octahedron-shaped bubble with six 
spins at its vertices fixed, while the spin sq in the middle is decoupled and plays no role any more. 



Obviously, the ratio ( |To| ) could be obtained by performing a usual canonical ensemble simulation of such a system 
with a defect, and measuring the probability distribution for its overall magnetization X)iefl5 Si ' ^ ne vame °f the 
ratio ([n]) would then be given by the ratio of the heights of the two neighboring bins in the corresponding histogram. 

In all cases of practical interest for the study of the scaling limit (sufficiently large systems, sufficiently small 
magnetization) the ratio ( JIOJ ) is close to 1. Otherwise a difference of one unit in the total magnetization would lead 
to a large change of probability: this would obviously be far from the scaling limit. Thus we always work with ft< 1. 

It is convenient to introduce a shorter notation Prs(x) = -Pc(X)iGflS Si = -Mc| s oi ' ' ' s 06) where x may be read as 
the magnetization of the remaining system if the factor N instead of TV — 7 makes a negligible difference, i.e. for large 
systems. We have to keep in mind that the notation Prs refers to a system with a defect whose type is not explicitly 
shown. Again restricting ourselves to local environments of the type 3+3—, we obtain 

(H) P+ P p Rs(M -jr)~ PnsjM + ±) 

{0) P+ + P- P RS (M-±) + P RS (M + ±y [ ] 



Thus we arrive at 



Also, due to h <C 1, Eq. (||) reduces to 



and we get 



/ \~ 1 1 dP Rs(M) 

vS ° } ~ NPrs(M) dM ■ {12 > 

(s ) = h, (13) 



-^logP^(M). (14) 



Defining the effective potential \Q S (M) (i.e., the Ginzburg - Landau fixed-M free energy) of the present system 
(with defect) by 

P R s(M) cx exp{-NV$ S) (M)}, (15) 

we get immediately 

- dV { ^ S) (M) , s 

h= _cS V 1 , 16 x 

dM v ' 

For large systems, the relative contribution of the defect is small, and thus h(M) is well approximated by V c g(M) for 
the finite system without a defect: 

P(M) cx exp{-JVF eff (M)}, (17) 

where the ellipsis stands for corrections vanishing at large N . 

As is well known, for finite 3D Ising models in a cubic box with periodic boundary conditions, the distribution 
P(M) has a double-peak structure MM at the critical point. Thus V a s (M) has a double- well shape, which immediately 
explains why h goes negative for small values of M. In Fig. [2] we show the quantitative comparison of h(M) (depicted 
by diamonds) and dV c s/dM (solid line). One observes that the correspondence between the points and the line clearly 
improves with increasing lattice size. To extract V a s from the Monte-Carlo-generated P(M) (Fig. |l|) we have exploited 
the fact, reported in fllO[], that for the system under consideration, P(M) can be well approximated by the ansatz 

f2 „ 2 / A#2 



P(M)aexp{-(^-l) (a M+ c)Y (19) 



which applies to the finite-size regime, i.e. the finite size is small compared to the bulk correlation length. We have 
fitted the Monte Carlo generated P(M) data accordingly determined the parameters a and c, and thus obtained 
dV e g(M)/dM in a simple polynomial form. 



It is also worth mentioning that the shape of P(M) for a given geometry (in our case, a cubic box with periodic 
boundaries) is universal at the critical point. That is, the parameters a and c have well-defined scaling limits when the 
system size grows to infinity. These values, a — 0.158(2), c = 0.776(2), have been determined in JlCJ ] by making use of 
a special model in the 3D Ising universality class, which has almost no corrections to scaling |ll| . The corresponding 
scaling form of dV e e(M)/dM is plotted by the dashed line in Fig. S. One observes that deviations from scaling 
(between the solid and the dashed lines) go down with increasing size, as they should. 

The results in Fig. g confirm the relation between the observable h(M) as defined above in the fixed-M ensemble, 
and the probability distribution P(M) in the canonical ensemble. The remaining discrepancy (between the diamonds 
and the solid line in Fig. g) is due to the "defect" discussed above. The question arises whether it is possible to 
modify our definition of h in order to suppress this discrepancy. We have found that this is indeed the case. Up to 
this point, we restricted ourselves to symmetric local environments (3+3—) to define h via Eq. (|3|). As has already 
been mentioned, using Eq. (H) one may use other types of local environments as well. In those cases the magnetization 
k = 2j=i s i enters the definition: 

(s )=t&nh(h + i3k). (20) 

Following the same arguments as before, we decompose the system in the local spin sq, its fixed neighbors, and the 
remaining system (RS) . This leads to the following generalization of Eq. (pd|) : 



, _ e? k P RS (M -£-£)- e-? k P RS {M - | + ^) 
^° ; e^P RS (M -&-£) + e-0*P RS (M £ + ±) 

fa tanh 3k * tt i ttt — — 

cosh 2 3kNp RS (M-&) dM 



A'/ 

~tanhfa-l ' dPR /W\ (21) 

V N P R s(M) dM J y ' 

Thus we arrive once again at Eqs. (OUTS). But we now have a different type of defect in the remaining system, and a 
shift of k/N in the magnetization of the remaining system; we neglect the latter effect. Now it seems plausible that one 
can suppress the influence of the defect by averaging over all configurations of the defect, weighted with their natural 
occurrence probabilities. Such an averaging should more faithfully reproduce the characteristics of a system without 
a defect. The modified determination of h is as follows. Sample configurations from the fixed-M ensemble. For each 
spin determine its orientation (+ or — ) and the type of its local environment (type 0, . . . , 6 for 0+6—, . . . , 6+0—). 
Accumulate these data by incrementing one out of 14 bins N q: +, N q —, where q = 0. . .6 denotes the type of local 
environment, and + or — denotes the local spin. The resulting population numbers satisfy X)o=o(-^9.+ + N q< -) = N. 
Then, for each q, find (s ) q = {N q ^ + — N q> _)/(N q<+ + A^_) and compute h q according to Eq. (||). Finally, 



improved = JZ ^X ' ( N 9,+ + N 1--)- ( 22 ) 

9 =0 



Applying this definition to our simulation data, we observe that, within the statistical accuracy, the discrepancy 
between h(M) and dV c g(M)/dM is indeed eliminated (Fig. g). 

IV. DISCUSSION AND CONCLUSIONS 

The relation (Il8|) looks exactly the same as the standard relation between the field and magnetization in the 
canonical ensemble: 

_ dV cS (M) 
11 ~ dM ■ [26) 

The observed differences between the properties of h(M) in the fixed-M ensemble and h(M) in the canonical ensemble, 
the most prominent of which is the nonmonotonic behavior of h(M) instead of the monotonic behavior of h(M), can 
be traced to the different definitions of the effective potential. The one that occurs in Eq. ( |i8[ ) is the fixed-M free 
energy, 



F eff (M) = -(1/N) \ogZ f (M), (24) 

while the one that enters Eq. (J23|) is defined via a Legendre transformation: 

V cS (M) = -(1/N) log Z c (h) + hM, (25) 

where 

M = (M) h (26) 

is the canonical average of the magnetization in an external field h. The partition functions Zf(M) and Z c (h) were 
defined in Section | Ina situation where fluctuations become negligible, the term hM in V e s cancels the field 
dependence of the Boltzmann weights. Then both definitions of the effective potential become equivalent, and both 
effective potentials approach the bulk form so that the difference between h and h vanishes. 

In a finite system, due to fluctuations, V c g differs from V e g. For instance, at the Ising critical point, the double-well 
form of V e s is absent: T4ff has a single-well form. 

Returning to Eq. (16), there is another finite-size effect: the difference between \Q S (M) and V e g(M) due to 



the presence of the defect. The relative contribution of the defect becomes small for large systems, and it is further 
suppressed by the improved definition of h, Eq. (p2|) . 

Thus for large systems and sufficiently high M , when the correlation length is small in comparison with the system 
size, the finite size effects are suppressed, the difference between V e g(M) and the bulk effective potential disappears, 
and our definition of h(M) reproduces the expected bulk behavior. 

In conclusion, we have studied the critical three-dimensional Ising model in the fixed-magnetization ensemble, in a 
cubic geometry with periodic boundary conditions. This was done by means of the recently developed geometric cluster 
Monte Carlo algorithm. We have defined a magnetic field-like observable h for this ensemble, studied its dependence 
on the magnetization M and explained its counter-intuitive nonmonotonic behavior: h first becomes negative and 
then positive with increasing M (Fig. ||). We have provided a quantitative description of h(M) by establishing a close 
relation with P(M) — the probability distribution of the magnetization in the canonical ensemble. The nonmonotonic 
behavior of h(M) can be understood as a manifestation of the same finite-size effect that is responsible for the double- 
peak shape of P(M) at the critical point. Furthermore we have shown that, when fluctuations are negligible, our 
definition reduces to the standard canonical relation M(h). Finally, we note that in the different context of the 
simulation of a system of particles whose number is fixed, a similar line of reasoning enables the determination of the 
chemical potential of the particles [12] . 
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FIG. 1. Top: Probability distributions P{M) of the magnetization per spin M — jj^2i i=1 s i> a ^ ^ ne critical point 
(/3 = 0.221654, h = 0), for the 3D Ising model, Eq. (hi), in a cubic box with periodic boundary conditions, for two lattice 
sizes: 12 3 (left) and 20 3 (right). The Monte Carlo data were obtained using the Swendsen-Wang cluster algorithm (720000 
configurations for each lattice size). The solid line is a fit according to Eq. (|l9|). For the 12 lattice, a — 0.268(13), c = 0.859(8), 
M = 0.3892(11). For the 20 3 lattice, a = 0.209(9), c = 0.839(11), M = 0.2984(9). Bottom: the difference between the data 
and the fit. 
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FIG. 2. Magnetic field h, computed as an observable in the fixed- M ensemble Eq. (0). The temperature, boundary conditions 
and lattice sizes are the same as in Fig. hi. Results obtained from Eq. (H), restricted to spins with a local environment of the 
type 3+3—, are shown as diamonds. Triangles correspond to the improved definition, Eq. (E3) . Solid lines show dV e g(M)/dM, 
where V e B (M) is exactly the same as in Fig. hi. The dashed lines show the universal shape of dV c s(M)/dM, using the universal 
(scaling-limit) values of the parameters: a = 0.158(2), c = 0.776(2) fllCf 



